{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecasting using Facebook's Prophet package\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "Forecasting is a fundamental data science task. Accurate forecasting is crucial when planning and allocating resources, setting goals and detecting anomalies. However, building reliable and robust forecasting models is challenging and often requires expert input. For example, there may be holidays or other 'one-off' effects that need to be taken in consideration for accurate forecasts. \n",
    "\n",
    "In this short tutorial we will provide an overview forecasting and discuss how to quickly implement your own forecast using the Prophet package recently released by Facebook Core Data Science team. The Prophet package fits an easily interpretable regression model which can account for seasonality trends (e.g., yearly or weekly effects) and holiday effects\n",
    "\n",
    "Throughout this tutorial we will be using an open access TFL dataset relating to the number of TFL cycles borrowed over a period of several years. This dataset will demonstrate considerable season effects (for example there are yearly seasonal effects as customers are likely to cycle in the summer) as well as and holiday effects (in particular, we will study the effects of bank holidays and the tube strike on the number of bikes borrowed!). This dataset will demonstrate considerable season effects (i.e., more likely to cycle in the summer) as well as trend and holiday effects. \n",
    "\n",
    "Before discussing the details of the prophet package, we provide an exploratory overview of the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "\n",
    "In order to get started, you will need to install the Prophet package. Simply use the following command on Mac and Linux:\n",
    "\n",
    "```\n",
    "pip install fbprophet\n",
    "```\n",
    "\n",
    "For Windows, follow the instructions from [this thread on the official development repository](https://github.com/facebook/prophet/issues/2#issuecomment-296321790)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import modules:\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from fbprophet import Prophet\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "#%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The dataset\n",
    "\n",
    "The dataset we will study corresponds to the total number of daily hires for Santander bikes over a period of approximately seven years. You can download the TFL cycle hire dataset from the following website: https://data.london.gov.uk/dataset/number-bicycle-hires"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in dataset:\n",
    "# Modify the directory name to the one you placed the dataset in\n",
    "#os.chdir('/somewhere/there/or/here/Data')\n",
    "dat = pd.read_excel('data/tfl-daily-cycle-hires.xls', sheetname='Data')\n",
    "\n",
    "# we are only interested in the first two columns\n",
    "dat = dat[ [u'Day', u'Number of Bicycle Hires']]\n",
    "dat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot( dat['Day'], dat['Number of Bicycle Hires'])\n",
    "plt.xlabel('Date')\n",
    "plt.ylabel('Number of cycle hires')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the figure above (of daily bike hires) we see several clear effects which are often present in forecasting challenges:\n",
    "\n",
    "1. Seasonality: There is a drop during the winter months and a rise during the summer months. This seasonality is easy to understand but needs to be accounted for in order to obtain reliable forecasts.\n",
    "2. Trends: There is a slight, overall upward trend. This means that the number of cycles hired has increased over time. \n",
    "3. Outliers: There are obvious outliers, notably during the summer of 2015. This coincides with the tube strike, which foreced commuters to find alternative routes. The histogram plot, shown below, highlights the presence of these large positive outliers. \n",
    " \n",
    "In the remainder of this tutorial we will introduce the prophet package and demonstrate how it may be used to handle each of the challenges discussed above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Histrogram of the data to demonstrate there are large positive outliers.\n",
    "# It is important to note that due to the seasonality, there may also be large \n",
    "# negative outliers, but we cannot spot them in the histogram below as \n",
    "# they are within the range of the data. \n",
    "#\n",
    "# the data appears to be approximately Gaussian, with a few outliers (we discuss these below).\n",
    "# so don't need to transform the data (e.g., via log transform)\n",
    "plt.hist((dat['Number of Bicycle Hires']), bins=40)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Prophet model\n",
    "\n",
    "The idea behind the prophet package is to decompose time series data into the following three components:\n",
    "\n",
    "- trends: these are non-periodic and systematic trends in the data,\n",
    "- seasonal effects: these are modelled as daily or yearly periodicities in the data (optionally also hourly), and\n",
    "- holiday/one-off effects: these are effectively outliers.\n",
    " \n",
    "Each of these components contribute *additively* to the observed time series. In other words, the number of cycle hires on any give day is *the sum* of the trend component, the seasonal component and the one-off effects. As a result, the model can be mathematically written as follows:\n",
    "\n",
    "$$ y(t) = g(t) + s(t) + h(t) + \\epsilon_t $$\n",
    "\n",
    "where $y(t)$ is the number of cycles hired on day $t$ and $g(t), s(t)$ and $h(t)$ correspond to the growth trend, seasonal effects and holiday effects respectively. Finally, $\\epsilon_t$ accounts for noise. \n",
    "\n",
    "We note that the approach taken in the Prophet packge is different from traditional forecasting models such as [ARMA](https://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model) models. For example, note that $y(t)$ is not a function of $y(t-1)$ but instead a function of the growth, season and holiday effects. The main advantages of taking an additive approach are as follows.\n",
    "\n",
    "- It provides us with more flexibility. For example, we can easily accomodate irregular measurements.\n",
    "- We can directly interpret each of the additive components and thus untangle the contributions of each aspect to the forecasts. This is important when we look to understand how the data was generated.\n",
    " \n",
    "We now discuss each of the components of the model\n",
    "\n",
    "### Trend component\n",
    "\n",
    "The default trend component in the Prophet package is a linear growth trend with some changepoints. If there are no changepoints, this means that the growth with continue to grow at some linear rate such that:\n",
    "\n",
    "$$ g(t) = \\alpha t$$\n",
    "\n",
    "Changepoints are introduced so that the linear rate of growth can vary. For example, it may be the case that there was 2% daily growth (corresponding to $g(t) = 1.02 \\times t$) for the first year and 1% daily growth (so that $g(t) = 1.01 \\times t$) afterwards.\n",
    "\n",
    "### Seasonal components\n",
    "\n",
    "Seasonality is a crucial component of many time-series and is clearly present in the TFL data we study here. The approach taken in the Prophet package is to consider various seasonal components with varying periodicities. For example, we may have one component with a periodicity of 7 days to capture weekly effects, and another with a periodicity of 365 days to capture yearly effects. By default, the Prophet package includes weekly and yearly seasonalities. This is what makes the most sense in this example, but it is also possible to include other seasonalities.\n",
    "\n",
    "### Holidays/outliers \n",
    "\n",
    "Finally, a well-document difficulty with time-series and forecasting data is the presence of outliers. These are data points which deviate significantly from the distribution of the remainder of the data. If these points are ignored they may possibly lead to poor forecasts and so it is important to consider these points carefully. Recall that in the case of the TFL cycle hire data, there were clear outliers around the time of the 2015 strikes (around the summer of 2015). Other important examples may include holidays or snow days.\n",
    "\n",
    "The Prophet package directly allows for such events by allowing the user to input a list of days which were holidays or potential outliers. The holiday component to the model is then:\n",
    "\n",
    "$$\n",
    "h(t) =\n",
    "\\begin{cases}\n",
    "\\kappa & \\mbox{if $t$ is a holiday/outlier}\\\\\n",
    "0 & \\mbox{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "In this way the scalar $\\kappa$ accounts for large positive/negative values of holidays or potential outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the model\n",
    "Now that we have the data ready as well as a high-level understanding of the forecasting model we will use, we are ready to fit our model!\n",
    "\n",
    "Before we begin, we must make sure we prepare the data in a correct format. We also need to prepare a `DataFrame` with holidays and the dates of the tube strikes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The prophet package expects input as a dataframe with the first column indicating time and \n",
    "# the second indicating the time series we wish to forecast\n",
    "\n",
    "dat['Day'] = pd.DatetimeIndex( dat['Day'] )\n",
    "\n",
    "# It also expects these columns to have the names 'ds' and 'y', so we rename them accordingly\n",
    "dat = dat.rename(columns={'Day': 'ds', 'Number of Bicycle Hires': 'y'}) \n",
    "dat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Following our discussion, we add a set of outliers/holidays for our Prophet model.\n",
    "# we make a separate DataFrame for bank holidays and for tube strikes\n",
    "\n",
    "# We get the list of bank holidays from the following csv file:\n",
    "\n",
    "bank_holidays = pd.DataFrame({\n",
    "    'holiday': 'BankHoliday',\n",
    "    'ds'     : pd.to_datetime( list(pd.read_csv('data/BankHolidayLists.csv')['Date']) )\n",
    "})\n",
    "\n",
    "strike_days = pd.DataFrame({\n",
    "    'holiday': 'strike',\n",
    "    'ds'     : pd.to_datetime( ['2017-08-05', '2017-08-06', '2017-02-06', '2015-07-09', '2015-07-08', '2015-03-08'] )\n",
    "})\n",
    "\n",
    "\n",
    "all_holidays_strikes = pd.concat( (bank_holidays, strike_days) )\n",
    "all_holidays_strikes.tail(n=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we are ready to fit a forecast model with prophet\n",
    "forecast_model = Prophet( growth='linear',  weekly_seasonality=3, yearly_seasonality=3, holidays=all_holidays_strikes )\n",
    "forecast_model.fit( dat )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting and understanding the model\n",
    "\n",
    "Now that we have build our forecasting model, we can interogate the model to understand what it is doing. \n",
    "The first way to do this is to see how the model fits existing data and what a forecast over 1 year looks like. This is shown below. From this plot we notice the model has done a good job of picking up the yearly seasonality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can now study the fit of the model - in order to do so, we need to creat another df\n",
    "df_dates = forecast_model.make_future_dataframe(periods=365, \n",
    "                                                include_history=True)\n",
    "model_predictions = forecast_model.predict( df_dates )\n",
    "plot_pred = forecast_model.plot( model_predictions )\n",
    "plt.legend(loc='best', fontsize=20)\n",
    "plot_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">\n",
    "In the plot above, the black dots correspond to the observed number of cycles hired each day (you'll notice that we don't have data beyond the beginning of 2018!). The dark blue line (labeled 'yhat') corresponds to the estimated number of cycles based on the estimated model. We notice it does a good job of capturing the main sources of variability in the data. Finally, the light blue lines correspond to the 80% confidence interval for the models predictions.\n",
    "</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get an even further understanding of the model, we can plot each of the model components. This is shown below. \n",
    "\n",
    "In the top panel we see the linear growth term. Recall that this term contained changepoints so that the rate of growth was allowed to vary over time. We notice that there is a large positive trend from around late 2010 until mid 2012. This is the initial adoption period. Then we notice a slight drop during 2013 before the trend stabilises. \n",
    "\n",
    "The second panel shows the effect of holidays and and tube strikes. The holidays have negative effects whilst the tube strikes have a positive effects. This suggests that the bikes are used for commuting.\n",
    "\n",
    "The final two panels show the estimated yearly and weekly trends of the model. As expected there is a large increase during the summer months when the weather is more pleasant. In the case of the weekly trend, we notice that there is a large drop during weekends. This also suggests that bikes are used for commuting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also study the model components - this allows us to further understand the data we study\n",
    "#\n",
    "# we notice the following:\n",
    "# initial trend during 2011-2013 which levels off (there appears to be a dip in summer of 2013, its unclear why)\n",
    "# clear yearly seasonlity with a peak in the summer and drop in winter\n",
    "# also a clear weekly seasonality - it seems the bikes are used mostly during the week (presumably for commuting) than on weekends\n",
    "# Finally, we note that bank holidays lead to a fall in the usage of the bikes, while the days of the strikes lead to significant increases!\n",
    "\n",
    "forecast_model.plot_components( forecast_model.predict(df_dates), \n",
    "                               uncertainty=False )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"color:blue\">\n",
    "**Note**: Finally, you may have noticed that the weekly plot appears as if it is continuous valued when we would expect it to be discrete (i.e., one value per day). Intuitively, it makes sense to fit the weekly component as an individual parameter per day of the week. This would result in a discrete weekly component, but it also requires us to estimate 7 parameters (one for each day of the week). This might not be the best way to procede. </span>\n",
    "\n",
    "<span style=\"color:blue\">\n",
    "In the case of the Prophet package, a slightly different approach is employed. Instead of estimate a parameter for each day as discussed above, the instead employ a [fourier series](https://en.wikipedia.org/wiki/Fourier_series) to model the daily effects. This effectively uses a periodic function to predict the number of cycles hired where the period is carefully adjusted to 7 days. For further details, see Section 3.2 of the [Prophet paper](http://amstat.tandfonline.com/doi/abs/10.1080/00031305.2017.1380080#.WoFjRpPFLdQ). \n",
    "</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations on fitting your forecasting model using the Prophet package! If you are interested in learning more about the underlying models and how they are estimated, you can find many more resources on the official Prophet website:\n",
    "\n",
    "https://research.fb.com/prophet-forecasting-at-scale/\n",
    "\n",
    "and the official starting guide on Github:\n",
    "\n",
    "https://facebook.github.io/prophet/docs/quick_start.html#python-api\n",
    "\n",
    "If you more interested in the model details and estimation, you can read the paper associated with the package here:\n",
    "\n",
    "https://peerj.com/preprints/3190.pdf"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
